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We formulate a real-space renormalization group (RG) approach for efficient numerical analysis 
of the low-temperature hopping dynamics in energy-disordered lattices. The approach explicitly 
relies on the time-scale separation of the trapping/escape dynamics. This time-scale separation 
allows to treat the hopping dynamics as a hierarchal process, an RG step being a transformation 
between the levels of the hierarchy. We apply the proposed RG approach to analyze hopping 
dynamics in one- and two-dimensional lattices with varying degree of energy disorder, and find 
the approach to be more accurate at low temperatures and computationally much faster, than the 
direct diagonalization. Applicability criteria of the proposed approach with respect to the time- 
scale separation and the maximum number of hierarchy levels are formulated. RG flows of energy 
distribution and preexponents of the Miller- Abrahams model are analyzed. 



I. INTRODUCTION 

Understanding the dynamics of diffusion and transport 
in systems/materials with disorder is of paramount im- 
portance in multiple branches of modern science and en- 
gineering. The areas where the disorder-related phenom- 
ena are ubiquitous include biology/ condensed matter 
physics, polymer and nano science. The nature of trans- 
port can be very different in various systems: excitons 
in molecular or semiconductor quantum dots aggregates 
and polymer chains/films; charge carriers (electrons or 
holes) in disordered semiconductors, phonons in disor- 
dered heat conductors. 

The nature of disorder can also vary a great deal from 
system to system. One of the most important sources 
of disorder is the energy dispersion - the dynamics of 
the system occurs in a configuration space with a rugged 
landscape, so that potentials barriers have to be over- 
come for e.g., an exciton to migrate within a polymer 
chain. Specifically, it is well documented that in a amor- 
phous conjugated polymers, excitons can be localized in 
short ordered segments of a polymer chain, where en- 
ergy of the exciton depends on the length of the segment 
giving rise to the energy disorder. An inevitable con- 
sequence is that the exciton, generated at an arbitrary 
site in the excitonic density of states distribution, will 
migrate (spectral diffusion) and, ultimately, tend to re- 
lax towards low-energy states.^ ^ It has to be noted that 
in amorphous films of conjugated polymers the energy 
disorder can be strong in a sense that the site-to-site 
energy variation can significantly exceed ksT at room 
temperature.^ ^ Possible mechanisms of exciton-excition 
coupling between polymer segments include Forster res- 
onance energy transfer (due to dipole-dipole Coulomb 
interaction),^ as well as Dexter mechanism due to the 
direct overlap of electronic wavefunctions.^ 

Transport in molecular/metallic nanocrys- 
tal/semiconductor quantum dot aggregates receives 
much attention lately because of multiple promising 
applications in photovoltaics, electronics and informa- 



tion processing. One of the most promising systems of 
this type is aggregates of semiconductor quantum dots, 
where the dispersion of energy of the lowest excited 
state comes from the quantum mechanical confinement 
of carriers in QDs with size distribution. In such a 
system, a photoinduced exciton has been demonstrated 
to efficiently migrate within an aggregatd^^^^. One of 
the possibilities to make quantum dots "talk" to each 
other, i.e., to make an exciton migrate, is via sol-gels - a 
special type of aggregates - which are of great pron iise 
for photovoltaic and thermoelectric applications.^^ 

Since the transport in the presence of disorder is typ- 
ically characterized by an absence of the well-defined 
time, size or energy scales, renormalization group (RG) 
is one of the most suitable methods to study dynam- 
ics under such conditions. RG has been extensively ap- 
plied be fore t o study the diffusion in energy-disordered 
systems ,^^^1^ however mostly for the case of the weak 
disorder. The case of the strong disorder, i.e., where the 
energy dispersion significantly exceeds the temperature 
(in energy units) of the environment, has been studied to 
a much lesser extent mostly due to the lack of appropri- 
ate methods. In this paper, we introduce the numerical 
RG method capable of studying the hopping diffusion in 
the presence of a strong energy disorder. 

The rest of the paper is organized as follows. The nota- 
tion and general theory is introduced in Section |TI[ The 
RG step (coarse-graining) is formulated in Section [ill 



General numerical results and the discussion of the RG 
flow are discussed in Sees. [IV| and [V| respectively. Sec- 
tion I VII concludes. 



II. HOPPING ON GRAPH 

We describe the hopping dynamics of a system with 
energy disorder as Markovian hopping between nodes 
of a graph, where the probability of hopping is defined 
by rate constants assigned to edges of the graph [see 
Fig. [TJa) for graph schematics]. Physically, each node 
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FIG. 1. (a) Schematic of a graph for hopping dynamics. Rate constants are shown for an edge connecting nodes 1 and 2. 
(b) Schematic of hopping dynamics on the graph with two local energy minima. At low temperature, the relaxation dynamics 
within each minimum is expected to be mush faster than the subsequent population transfer between the minima. 



of the graph can represent, for instance, a random- length 
ordered (e.g., straight) segment of a polymer chain with 
kinks, so that an exciton is localized within the segment. 
Accordingly, a population of a node represents a prob- 
ability to find an exciton at the corresponding segment. 
A scalar parameter, energy, is assigned to each node of 
the graph and two hopping rate constants (forward and 
backward) are associated with each edge of the graph, as 
illustrated in Fig. [TJa) for nodes 1 and 2. The important 
characteristic of the so formed graph is its connectivity, 
which, in the current context, can be defined as the num- 
ber of edges connected to a node. In general, a hopping 
rate constant can be non-zero between any two nodes 
of the graph, however distant they are. Physically, this 
corresponds to a situation where, e.g., a non- vanishing 
exciton-exciton coupling is present between two very dis- 
tant segments of a polymer chain, or between distant 
quantum dots within an aggregate. Practically, however, 
hopping rate constants typically decay rapidly with dis- 
tance, i?, e.g., as ~ Rq/R^ in the case of the dipole- 
dipole (Forster) or ~ ^-R/Ro \ji the case of the exchange- 
type (Dexter) interaction, where Rq is some characteris- 
tic range of interaction. Henceforth, we allow edges to 
link only nearby nodes, thus restricting the graph con- 
nectivity to be finite even in the thermodynamic limit. 

The hopping dynamics on the so introduced graph is 
governed by the following kinetic equation 



dt\p{t)) = w\p{t)), 



(1) 



where \p{t)) is the vector of time-dependent populations, 
i.e., Piit) = [\p{t))]i^ i = 1,...,A/' {N is the number of 
nodes in the graph) is the population of the node i at 
time t. The rate matrix is constructed as follows.— Its 
nondiagonal elements read as 



(2) 



where kij is the rate constant for hopping from node j 
to node i fsee Fig.jlja)]. Once the non-diagonal elements 
are defined, the diagonal ones are determined from the 
population conservation, which requires the sum of ele- 



ments in any column of W to vanish identically, i.e., 

Wii + J2Wji=0- (3) 

Therefore, Eq. ([T]) is the gain- loss equation for popula- 
tions of nodes. 

An additional constraint we impose on the rate con- 
stants is the detailed balance, i.e.. 



h ■ 



(4) 



where j3 = l/ksT = 1/T {ks = 1 throughout the paper). 
Energy of node i is denoted by Ei. Eq. Q guarantees 
that state [|Peg)]i ^ e~^^^ is the equilibrium state of the 
system, i.e., it is the steady state {dt\p{t)) = 0) with no 
currents^ 

Eq. can be formally solved as \p{t)) = e^^\po)^ 
where \po) = \p{t = 0)) is the vector of initial popula- 
tions. The rate matrix, which satisfies Eqs. (|3| and Q, 
can be diagonalized as W\n) = — A^ln),'^ yieldin^^ 



E 



n){n\po) 



(5) 



Because of the detailed balance, the lowest eigenvalue is 
always zero, corresponding to the equilibrium. 

Eq. ([5| completely solves the problem introduced by 
Eq. (Tj), since any observable which depends on node 
populations at an arbitrary moment of time can now be 
expressed via Eq. ([5|. For example, the average energy 
of the system at time t is given by 

^E{t)) = ^-^'P(^)^ = T.n{E\n){n\Po)e-'"\ 



(ibW) (ibo) 

where {E\ and (1| is vectors of node energies and units, 
respectively. Deriving Eq. ([6|, we explicitly used the pop- 
ulation conservation requirement. 

III. COARSE GRAINING 

Eq. (|5| gives the exact solution of the problem of hop- 
ping dynamics on graph. However, this solution might 



3 



not be optimal for the following reasons. Dispersion of 
node energies leads to the emergence of local minima 
where the population can be trapped at low tempera- 
ture. From the eigenproblem perspective, it results in 
a very broad eigenspectrum and coexistence of differ- 
ent timescales: the relaxation within local minima occurs 
rapidly, but the population transfer between local minima 
occurs slowly due to the necessity to overcome potential 
barriers. This is illustrated in Fig.jljb), where relaxation 
within local minima (1-2-3 and 5-6-7) can be fast, but it 
is hard to overcome the potential barrier (3-4-5), which 
can result in a very slow inter-minimum transfer at low 
temperatures. The natural criterion for the separation of 
timescales for these two processes can be introduced as 
7 = SE/T > 1, where SE = {E'^ - {E)Y^'^' At T ^ SE 
the eigenvalues of W can be spread over many orders of 
magnitude, so the numerical diagonalization of the rate 
matrix becomes inaccurate due to numerical round-off 
errors. Furthermore, the existence of the small param- 
eter, i.e., 1/7, might simplify the solution and analysis 
of the problem and/or allow for much larger systems to 
be dealt with than the brute-force direct numerical di- 
agonalization allows. Similarly to how the introduction 
of a small parameter simplifies a problem in quantum 
mechanics by allowing for a perturbative solution, we in- 
tend to explicitly use the smallness of I/7 to simplify 
the solution of the problem and make it more physically 
transparent. 

The straightforward idea is to explicitly separate the 
time scales by first solving for the population relaxation 
dynamics within each local minima. Then, the inter- 
minima transfer problem is solved considering population 
in each minima to be in its own quasi-equilibrium. This is 
similar to the transition state theory, where the timescale 
separation allows to treat both reactants and products 
as equilibrated within respective potential wells. Using 
the analogy between energy dispersion of the graph and 
rain/water collection in rugged landscape we will refer 
to local minima as to drainage basins or simply basins. 
Thus, we will split the entire graph into a collection of 
basins and, first, solve for the population dynamics within 
each basin independently. Then, the population is con- 
sidered to be equilibrated within each basin and dynam- 
ics of population transfer between basins will be solved 
for. To make a step further, one can notice that the 
population transfer between the equilibrated basins can 
be again treated as a Markovian hopping process on a 
new graph, which describes inter-basin hopping dynam- 
ics. This new graph can be thought of as a result of 
coarse-graining of the initial lattice. The construction of 
this new graph out of the initial one can be thought of as 
a renormalization group (RG) step (or transformation). 
This RG transformation can be performed multiple times 
successively giving rise to the hierarchy of levels of coarse- 
graining. At each such level, the population is being 
transferred between the basins, each basin consisting of 
nodes of the graph belonging to the previous hierarchy 
level. Fig. [2] schematically depicts the consecutive RG 



steps applied to the original (microscopic) graph. In this 
figure, the original (e.g., microscopic) graph is designated 
by Go- The basins are outlined by the red contours, and 
upon the RG step, each basin within Gobecomes just a 
single node (red circle) in a new graph, Gi. In turn, 
the RG transformation can be applied to the new graph, 
yielding G2 etc. In what follows, each graph within the 
RG hierarchy will be designated by G^, where index h 
denotes the number of RG transformation applied con- 
secutively to the original graph Gq. The ordered set 
(lowest to highest) of the exact eigenvalues of graph Gh 
(e.g., obtained via numerical diagonalization) is denoted 
by = {Ai, A2, Aat}. The set of eigenvalues corre- 
sponding to relaxation within basins of Gh is denoted 
by A'^. Further, we introduce a notation for the entire 
spectrum of the system, obtained by H successive RG 
transformations, as 

H-l 

Sh= [j K^Ah. (7) 

h=0 

Then, if the RG transformation were exact then the fol- 
lowing equality would hold 

So = Si = S2 = (8) 

However, since an RG step is approximate, the validity 
of this expression depends on the magnitude of 7 at each 
hierarchy level. In the remainder of the section we discuss 
the specific procedure of coarse-graining. 

A. Single RG step 

The following algorithm of coarse-graining is used: 

1. In the entire graph, identify nodes which are linked 
only to nodes of higher energy. Each such low- 
energy node will be called a "seed" of a basin. 

2. For any other node (of index z), a node j is found 
so that the following requirements are satisfied: (i) 
there exists an edge linking nodes i and j, (ii) Ei > 
Ej and (iii) of all indices j satisfying the first two 
requirements we choose the one maximizing kji. If 
node j has been associated with any basin (e.g., 
node j is a seed), we associate node i with the same 
basin. 

3. Repeat Step 2 until all the the nodes of the graph 
are associated with basins. 

The first execution of Step 2 will associate with basins 
only the nodes connected directly to basin seeds. Then, 
nodes connected to seeds through a single node (i.e., via 
two edges) will be associated with basins etc. At very low 
temperature (i.e., 7 1), this approach guarantees that 
a system, initially residing at a particular node, will pref- 
erentially relax to the seed of the basin the initial node 
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belongs to. Here, "preferential" means that we construct 
the basins using maximum kij requirement, and, there- 
fore, we neglect e.g., the second largest kij which can 
lead to non-negligible (but typically small) relaxation of 
the exciton to other seeds. This can be illustrated by 
initially putting an exciton to node 4 of the graph shown 
in Fig. [TJb). As seen, it can relax to both seeds 2 and 6, 
but we associate node 4 with the left basin if ks4 > /C54, 
and with the right basin if ks4 < /C54. 

Once basins are formed, we construct rate matrices 
(b is the basin index) for each basin as if it is closed, i.e., 
we set kij to zero if nodes i and j belong to different 
basins. Then, a typically small rate matrix of each basin 
is is numerically diagonalized 



W''\n,b) = -\l\n,b), 



(9) 



where the lowest eigenvalue vanishes exactly (i.e., Aq = 0) 
because of the detailed balance, just like in the case of 
the entire graph (see above). The eigenvectors, corre- 
sponding to these quasi-equilibrium eigenmodes (one per 
basin), read as (right and left eigenvector, respectively) 



ieb 



K0,6|], = [(l6|]i=5ie6, 

where ^^^5 equals 1 if node i belongs to basin h and zero 
otherwise. Partition function of basin h is given by Z5 = 



'}Zieb^~^^^ ' straightforward to check that these 

eigenvectors are orthonormal, i.e., (0,6'|0,6) = (^55/. 

The problem of "preferential" relaxation described 
above can lead to a significant underestimation of rate 
constants corresponding to intra-basin relaxation, 
(n 7^ 0). Indeed, at very low temperatures a popula- 
tion of node 4 (Fig. [T]) would decay with rate ks4 + /C54. 
Associating node 4 with a specific basin (left or right), 
as it is done in the coarse-graining procedure described 
above, would allow it to decay with rate max{ks4^ks4) 



which can introduce a significant error, up to factor of C 
- graph connectivity. To correct for this, we recalculate 
eigenvalues of intra-basin rate constants via 



A^ 



-{n,b\W\n,b). 



(10) 



Using the language of perturbation theory in quantum 
mechanics, on can say these new eigenvalues are first- 
order corrected relative to the zeroth-order approxima- 
tion given by Eq. Eigenvalues A^ (n / 0) constitute 
A'^ at a given hierarchy level /i, and are not modified any 
more by any further RG transformations. 

The rate matrix for the the coarse-grained graph, de- 
scribing the population transfer between basins, is con- 
structed as follows. The rate constant for population 
transfer between equilibrated basins is defined as 



kb,i, = {o,b'\w\o,b) = Y,Y.'' 

ieb' jeb 



(11) 



Physically, this means that we take equilibrium popula- 
tions of basin h and multiply them by rates of populations 
transfer to basin h' . It is easy to see that 



h'b/h 



bb' 



-/3(F,,-F5) 



(12) 



where the free energy of a basin is defined as F5 = 
— TlnZ^. Thus, the new graph - the result of RG trans- 
formation applied to the graph of the previous hierarchy 
level - is constructed as: 

1. Each node is the basin of the previous graph. 

2. Energy of each node is the free energy of the equi- 
librated basin. 

3. Rate constants of population transfer between 
nodes are given by Eq. (11) 



4. Initial populations of each node are given by 
{0,6|po). 
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5. Time-dependent population of each node is inter- 
preted as the total population of the equilibrated 
basin. 



B. RG representation in full dimensionality 

In contrast to standard RG procedures in statistical 
physics, where upon an RG step the microscopic informa- 
tion is discarded as irrelevant, we keep the entire eigen- 
spectrum during the coarse-graining procedure. In the 
coarse-graining procedure we proposed, the rate matrix 
W is partially diagonalized at each RG step in a sense 
that the first RG step applied to W transforms the rate 
matrix to a block-diagonal one, where one of these blocks 
is exactly diagonal itself, as it consists of eigenvaluesA^ 
(i ^ 0), corresponding to the intra-basin relaxation. The 
other block, non-diagonal within, is the rate matrix for 
inter-basin transitions, constructed out of rates defined 
by Eq. (11). Each successive RG step further decimate 



(MA) modeP 



this internally non-diagonal block and correspondingly 
increases the size of the internally diagonal one keeping 
the overall dimensionality constant. In other words, the 
RG approach introduced here is a method to approxi- 
mately diagonalize a certain class of matrices - rate ma- 
trix W at low temperature. It is, therefore, clear that 
Eq. ([5| is still valid as long as the RG-based diagonal- 
ization of W is accurate. For example, Eq. ^ can still 
be used to find the time-resolved average energy of the 
system. 

Clearly, only finite number of successive RG transfor- 
mations can be applied to a finite system, since at each 
successive level of the hierarchy, the size of the graph 
is rapidly decreasing (see below). In fact, after certain 
number of successive RG step the direct diagonalization 
it typically possible even for a very large initial system. 
The validity and accuracy of coarse-graining procedure 
can impose an additional constraint on the maximum 
number of RG steps. Indeed, 7 may become not too large 
at a certain level of hierarchy effectively rendering the RG 
step inaccurate (see detailed discussion in Sec. M). 



IV. NUMERICAL RESULTS FOR 
COARSE-GRAINING 

To demonstrate the validity of the proposed RG ap- 
proach at low temperatures, we introduce the following 
generic system with energy disorder. The regular lattice, 
one-dimensional chain (ID) or two-dimensional square 
lattice (2D), is constructed with a random energy as- 
signed to each site via sampling the standard normal 
distribution {{E) = 0, SE = {E'^) = 1). Only near- 
est neighbors are linked (with periodic boundary condi- 
tions assumed) , thus forming graph Gq with connectivity 
(7 = 2 and C = 4 in the ID and 2D case, respectively. 
Hopping rate constants are assigned to each edge of the 
graph using the modified version of the Miller- Abrahams 



J /cO. exp(-/3A^,,), AE,j > (up - hop). 



AEij < (down - hop). 



Here, a single preexponent is assigned to each edge 
k^j = k^j^ = exp(— /3e^), 



(13) 



(14) 



where e^j is an activation energy of the hopping along 
edge ij. Activation energy is assigned via sampling a 
uniform distribution so that e'^j G (0, 1). The presence of 



this temperature-dependent preexponent makes Eq. (13) 
different from the classical MA model, the latter recov- 
ered by fixing the activation energy, e.g., efj = 0. The na- 
ture and implications of this modifications are explained 
below. In the rest of the section numerical results are 
presented for ID and 2D systems with rates specified by 
Eq. ( pl) 

Fig. [3] presents the results RG-based calculations of 
eigenvalues for a ID system with 5000 nodes. Panel 
(a) shows the eigenspectra for the the direct diagonal- 
ization (Ao) and the RG approach with a coarse-graining 
(Aq, Ai) at a temperature equal to the magnitude of 
the energy disorder (T = SE = 1). As expected, at 
7 = T/SE > 1 the performance of the RG approach is 
poor. Specifically, it underestimates the large eigenval- 
ues (i.e., intra-basin relaxation rates) since an exciton 
can escape a basin (and therefore increase an eigenvalue 
of a transient mode) on the same timescale as the intra- 
basin relaxation occurs. In turn, not accounting for such 
basin escape processes at short times (Aq) leads to their 
emergence at larger timescales (Ai), resulting in overes- 
timations of smaller eigenvalues. 

Numerical results for the significantly smaller temper- 
ature, T = 0.05, are shown in panels (b) and (c). As 
seen, the performance of the RG method has been im- 
proved, as compared to the large temperature case, as 
it nearly coincides with the exact diagonalization results 
for a large (higher) part of the spectrum. In fact, at such 
low temperatures the RG approach outperforms the di- 
rect diagonalization. The kink seen in the DD spectrum 
at ~ 10~^^, followed by a rapid fall, is due to numerical 
round-off errors: ~ 220 eigenvalues turned out to be neg- 
ative in Aq. In comparison, a single coarse-graining gives 
only ^ 150 negative eigenvalues in Ai [see panel (b)], 
which is seen by a lower position of the kink. Further- 
more, two successive RG steps give only ^ 50 negative 
eigenvalues in A2, as is seen in panel (c) by the absence 
of the kink and very low onset of the rapid fall for the 
RG results. 

Numerical results for the 2D square lattice with 50 x 
50 = 2500 nodes are shown in Fig. [4j Similarly to the ID 
case (Fig. |3|, the RG approach performs poorly when 7 
is not large, as is shown in panel (a). Panel (b) demon- 
strates the much better performance at lower tempera- 
ture, where the RG transformation has been applied four 
consecutive times. The numerical round-off error is again 
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FIG. 3. Comparison of the RG approach with the direct diagonahzation (DD) for the ID chain, Aq. (a) DD versus a single 
RG step at high temperature (T = 1). (b) DD versus a single RG step at low temperature (T = 0.05). (c) DD versus two RG 
steps at low temperature (T = 0.05). 





FIG. 4. Comparison of the RG approach and DD for the 2D square lattice of size 50 x 50 = 2500. (a) DD versus a single RG 
step at T = 1. (b) DD versus 4 RG steps at T = 0.05. (c) CPU time versus number of RG steps. 



seen as a vertical drop of the DD spectrum in the low part 
of the spectrum. One can also notice that at the same 
low temperature and the number of nodes in the lattice, 
the ID system [Fig. [sj^c)] possesses a somewhat broader 
spectrum than the 2D one [FigQb)]. The reason for 
this is that the ID diffusion is an exact lower boundary 
for nD diffusion, since hopping between any two nodes 
in nD can always be considered as hopping in ID but 
with the "shortest" path, i.e., the path with highest rate 
constants. 



Fig. ^c) shows the computational time required to ob- 
tain the entire spectrum of hopping dynamics with H 
RG steps. The hierarchal RG procedure has been imple- 
mented into Fortran 90 codes, and all the calculations has 
been performed at Intel(R) Xeon(R) CPU (2.66GHz). As 
is seen, the DD method is not only plagued by the nu- 
merical round-off errors at low temperatures, but also 
by far the slowest one, since one large matrix has to be 
diagonalized at once. In contrast, at each step of RG 
many small matrices has to be diagonalized, which enor- 
mously speeds up calculations. Furthermore, since di- 
agonahzation of small matrices (intra-basin relaxation) 
is absolutely independent, the efficient parallelization is 
possible. 



Time-resolved luminescence 



To illustrate how the proposed RG approach can be 
used to evaluate various spectroscopic observables, we 
simulate the process of exciton migration within the 
energy-disordered lattice and calculate the time-resolved 
luminescence corresponding to this process. Specifically, 
we evaluate the time-resolved average luminescence en- 
ergy. If an exciton is localized at node i (of energy Ei), 
it can recombine radiatively emitting a photon with en- 
ergy = Ei. Therefore, if we assume for simplicity 
that the rate of radiative recombination does not depend 
on exciton energy, then the average luminescence energy 
at time t is given by {huj){t) = {E\p{t)), i.e., by Eq. (|6| 
for the average energy of the system at time t. As uni- 
versally in this work, we assume SE = 1 and {E) = 0, 
the latter equality implying that {E\p{t)) we calculate is 
not the absolute average of an exciton energy, but in- 
stead the time-dependent average energy measured rel- 
ative to the average node energy of the system. The 
time-dependent average energy of the 2D system of 2500 
sites is shown in Fig. |5] For definiteness, all the states 
are assumed to be populated with equal probability at 
t = 0, i.e., \p{t = 0)) = This assumption 

results in {E{t = 0)) = for both the exact solution 
obtained via direct diagonahzation. So, as well as ap- 
proximate RG-based ones, Sh {h > 0). At short times, 
the agreement between exact and approximate solutions 
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FIG. 5. Average energy of the graph versus time for the graph of size 50^ = 2500. Temperature is set to T = 1 and T — 0.2 
in panels (a) and (b), respectively. 



is still good in panel (a), but at longer times {t ^ 1 — 10) 
the RG approach is inaccurate, as expected for this high- 
temperature case, T = 1. Within the same time domain, 
the agreement of the RG approach and the direct diago- 
nalization is better in panel (b) , where temperature is sig- 
nificantly lower (T = 0.2), and, therefore, RG is expected 
to produce accurate results. At large times the system 
is equilibrated in panel (a), so RG and the direct diago- 
nalization results coincide again. At lower temperatures, 
panel (b), the equilibration of system is taking signifi- 
cantly longer, and has not yet been reached at maximum 
observation time of t = 10^. 




V. RENORMALIZATION GROUP FLOW 

The analysis of dynamics of a system on different time- 
and length-scales, as well as interaction between these 
scales, can often be conveniently done by studying the 
RG flow - change of parameters of the system (i.e., hop- 
ping graph in this work) under successive RG transforma- 
tions. This can be illustrated by discussing the RG flow 
of the energy distribution of graph nodes. Fig. [6] shows 
the distribution of the original microscopic graph, Go, 
which is obviously just the standard normal distribution 
{{E) = 0, SE = 1), as well as the energy distributions for 
graphs obtained from Go by RG transformations: Gi, 
G2 and G3. As is seen, the node energy distribution 
flows to lower mean energies and lesser energy disper- 
sion, i.e., the distribution grows sharper under successive 
RG transformations. The extracted mean energies and 
standard deviations are plotted by black circles and red 
squares, respectively, in the inset. The reduction of the 
mean energy can be easily understood from the fact that 
at very low temperatures, energies of nodes of the new 
graph are the energies of basin seeds, which are found by 
requiring each such point to be energetically lower than 
its immediate neighbors. The reduction of the energy 
dispersion upon the RG step is less trivial and can be 
understood from order statistics. Speciflcally, the mean 
and the standard deviation of a minimum of n sample 
points (in the limit of large n), each obtained by sam- 



FIG. 6. Energy distribution of nodes corresponding to hier- 
archy levels from h — (initial graph) to h — 3 for the 2D 
system consisting of 200 x 200 = 40000 nodes at T = 0.05. 
The inset shows mean (black circles) and the standard devi- 
ation (red squares) corresponding to the distributions in the 
main figures. The size of the graphs for the each hierarchy 
level are 40000, 8046, 1123 and 142 for Go, Gi, G2 and G3, 
respectively. 



pling the standard normal distribution is approximately 
given hy '^'^ 



{E): 

5E'^ 



(21ogn)-^/^ 



(15) 



At low temperature, node energies of a graph obtained 
by the RG transformation are equal to the energies of 
the seed nodes of the preceding graph, which are, in 
turn, obtained by essentially an order statistics - out of 
each immediate neighborhood (of size n) the node of the 
minimum energy was identifled. These expressions en- 
code the reduction of both mean and dispersion of the 
node energy distribution upon an RG step, which is in 
qualitative agreement with the numerical results (inset 
of Fig. [6|. Quantitatively, the average basin size in the 
original graph is n = ^ 5, so that Eq. (15) yields 

5E ^0.56 and (E) ^ —1.8. The numerical results (from 
the inset) are SE ^ 0.67 and (E) ^ —1.2, so the agree- 
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weak RG 
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V 



7 = 1 



Diffusive fixed point 

(7 = 0) 



Infinite disorder 

(7 = 00) 



FIG. 7. Schematic representation of the RG flow of 7 
5E/T. 



ment is rather good considering the asymptotic (large n) 
nature of Eq. (15). 

Fig. [t] schematically depicts the RG flow of 7 = 5E/T^ 
i.e., how 7 changes when successive RG transformations 



are applied. The RG approach proposed in this work is 
valid only at 7 ^ 1, as is shown earlier, and, therefore, 
one can call it strong disorder RG. Fig. |6] and Eq. (15) 
suggest that upon successive RG transformations, 7 flows 
to lower magnitudes, as illustrated in Fig. [7| This im- 
poses a natural restriction on the maximum number of 
valid RG transformation, since above certain number of 
transformations 7 becomes ~ 1, and the strong disorder 
RG becomes invalid. 

On the other hand, the weak disorder RG, discussed 
by Deem and Chandlei^^ demonstrates the flow in the 
same direction, i.e., reduction of 7, ultimately reaching 
the stable fixed point at 7 = 0, which is the pure dif- 
fusion. Neither of these RG approaches is valid in the 
intermediate disorder regime, where 7^1. However, the 
exact solution for the ID model suggests that the long- 
time behavior of ID model at arb itrari ly strong gaussian 
disorder is always purely diffusive j ^ ^ I ^ ^ I Fur t her more , this 
exact solution for the ID model is an exact lower bound- 
ary for models with higher dimensionality.^^ Therefore, 
observing hopping in a model with gaussian energy dis- 
order at increasing lenghtscale is expected to always cor- 
respond to a decreasing effective disorder. We therefore 
put dotted arrow in Fig. [7| to reflect the expected flow 
direction at 7 ~ 1, even though we are not aware of any 
RG approach valid in this region. 

Finally, we briefly discuss the RG flow of the distri- 
bution of the preexponents of the modified MA model. 
Within the standard MA model preexponents are all de- 
generate. It is however quite clear that even with this 
degeneracy present, the preexponents of the rate con- 
stants obtained upon the first RG step will acquire cer- 
tain finite dispersion since preexponents of inter-basin 
rate constants are determined by activation barriers be- 
tween basins, and the heights of these activation barri- 
ers is a random variable. Thus, it is expected that the 
degeneracy of the preexponents is irrelevant in the RG 
sense, i.e., the system flows away from the preexponent 
degeneracy upon successive RG transformations. To il- 
lustrate this, we compare the eigenspectra of the modi- 
fied and standard MA models, the latter obtained from 
the modified one by setting the activation energies to 
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FIG. 8. The eigenspectrum for the modified MA model 
[tij ^(0,1), black line), as well as for the standard MA 
with the fixed preexponent obtained from Eq. ( 14 ) by fixing 
the activation energy: e^j — (blue squares), e^j — 0.6 (red 
circles) and e^^ = 1 (magenta diamonds). All the eigenspectra 
are evaluated for the 2D system with 2500 nodes at T = 0.05 
applying two successive RG steps. 



a constant value. Fig. [s] shows these eigenspectra {S2 at 
T = 0.05 for the 2D system with 2500 nodes) for the mod- 
ified MA model (black line) with e^^ sampled uniformly 
within (0, 1) interval, as well as well as eigenspectra for 
the standard MA with activation energies set to to e^^ = 1 
(magenta diamonds), e^^- = 0.6 (red circles) and e^^- = 
(blue squares). One can see that in the high- A portion 
of the eigenspectrum modified and standard MA models 
produce qualitatively different results. Specifically, the 
standard MA model produces almost flat eigenspectrum 
due to the high degree of degeneracy of eigenvalues. How- 
ever, the lower part of the spectrum is qualitatively the 
same, except for a vertical shift. The origin of this shift 
is clear, e.g., at e^^- = 1 all the rates are lower than in the 
case where activation energy is distributed between and 
1, and, therefore, the spectrum with this fixed activation 
energy is expected to be generally lower than the one with 
finite dispersion of activation energies (compare the black 
line and magenta diamonds in Fig.|8|. Similar argument 
explains why the eigenspectrum with e^j = is generally 
higher than the one with the dispersion (compare the 
black line and blue squares in Fig. [8| Furthermore, it is 
possible to numerically find the fixed activation energy, 
e'^j =0.6, to make the both models practically coincide in 
the lower part of the spectrum. Therefore, the degener- 
acy of preexponents is irrelevant in the RG sense, which 
is reason why we adopted the modified MA model from 
the onset. 



VI. CONCLUSION 

In this paper, we formulate and computationally imple- 
ment a numerical RG approach to study the hopping dy- 
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namics in energy- disordered lattices at low temperature. 
Our approach is somewhat similar to the so called "strong 
disorder RC'^^M^, which has been developed originally 
for disordered Ising models. This similarity stems from 
the fact that in both cases fast degrees of freedom can 
be treated independently from the slow ones, giving rise 
to renormalization of the latter. We emphasize however, 
that in this work the fast degrees of freedom are not 
actually "integrated out", since the entire spectrum is 
recovered, not just its low-energy portion. The other dif- 
ference is that the formulated RG approach is dynamical 
in a sense that the dynamics in the form of eigenvalues 
and eigenvectors is recovered at all timescales, and not 
just equilibrium properties of a system of interest. 

In this work, the formulated RG approach is applied 
to discrete systems, represented by graphs. In fact, in 
the case of strong disorder the approach can be straight- 
forwardly applied to continuous systems. Indeed, even if 
the microscopic system is initially continuous, the popu- 
lation trapping in local minima leads to the appearance 



of the discrete-like dynamics, i.e., hopping between now 
continuous basins. In other words, the RG transforma- 
tion can still be introduced, and its first application to 
the continuous system in the strong disorder regime im- 
mediately leads to discrete hopping on a graph. Thus, 
the continuity of a system is irrelevant in the RG sense 
at strong disorder. Interestingly, it is exactly the op- 
posite in the weak disorder limit, where lattice becomes 
irrelevant and continuity "sets it" finally leading to pure 
diffusion at large timescales. 

We expect the developed approach to become useful in 
studying hopping dynamics in various systems of prac- 
tical interest including polymer films and nanoparticle 
aggregates. 
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